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ABSTRACT 



The Naval Postgraduate School (NPS) has actively explored the design and imple- 
mentation of networked, real time, three-dimensional battlefield simulations on low 
cost, commercially available graphics workstations. The most recent system, NPSNET, 
has improved in functionality to such an extent, that it is considered a low cost version 
of the Defense Advanced Research Project Agency’s (DARPA) SIMNET system. In 
order to reach that level, it was necessary to economize in certain areas of the code so 
that real time performance occurred at an acceptable level. One of those areas was in 
aircraft dynamics. However, with "off-the-shelf computers becoming faster and 
cheaper, real-time and realistic dynamics ate no longer an expensive option. The real- 
istic behavior can now be enhanced through the incorporation of an aerodynamic mod- 
el. To accomplish this task, a prototype flight simulator was built that is capable of 
simulating numerous types of aircraft simultaneously within a virtual world. Beside be- 
ing easily incorporated into NPSNET, such a simulator will also provide the base func- 
tionality for the creation of a general purpose aerodynamic simulator that is particularly 
useful to aerodynamic students for graphically analyzing differing aircraft's stability 
and control characteristics. This system is designed for use on a Silicon Graphics work- 
station and uses the GL libraries. 
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I. INTRODUCTION 



A. MOTIVATION 

The current state of the an in simulation technology has provided today's military with 
many extremely valuable training experiences that could not have been obtained elsewhere 
and, as a result, has greatly increased survivability and readiness. From flight simulators, 
which allow a pilot to explore the edge of the flight envelope without endangering crew or 
multi-million dollar assets, to battlefield simulators, which allow' entire fighting divisions 
to practice command and control without having to incur the enonnous costs of running a 
full blown field exercise, computer simulation has become a vay of doing business within 
the military. 

One simulation system designed by the Defense Advanced Research Projects Agency 
(DARPA) is the Simulation Networking (SIMNET) [Thorpe, 1987). S1MNET is a 
networked battlefield simulator that allows multiple use interaction on the battlefield at 
many different levels. Vehicle simulators, such as tanks and aircraft, connect to the network 
and become part of a three dimensional world. On the lowest level, this system allows the 
vehicle drivers a chance to study defensive and offensive tactics. On higher levels, 
Company, Battalion and Regiment Commanders can experiment with the movement of 
forces, thereby learning what will and won’t work at their particular level. 

At the Naval Postgraduate School (NPS), an effort to develop a SIMNET type system 
based on commercially available, general purpose, graphic workstations, has been active 
for a number of years. This system, NPSNET, consists of Silicon Graphics workstations 
attached to a local area Ethernet [Zyda.et al.,1992). Eventually, NPSNET will become a 
node on the SIMNET network. 
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Over the past few years, as workstations became available that were both faster and 
cheaper, capabilities, that were once too computationally expensive to incorporate, have 
become feasible. One such feature is realistic vehicle dynamics. As a result it is now 
possible to enhance the behavior of the system by the incorporation of realistic vehicle 
dynamics. A ground vehicle dynamics model currently in development at NPS will provide 
a greatly improved system for ground vehicle interaction with the enviomment. Because 
aircraft dynamics differs in many respects to ground vehicle dynamics, a separate dynamic 
and orientation model for air vehicles is required. This work provides such a model for 
incorporation into NPSNET. 

Most of the research and development of aerodynamic simulators has involved either 
the development of large multi-million single aircraft platforms utilizing specialized 
computer hardware to increase processing speed, or the development highly proprietary 
flight simulation programs for use on a personal computers. As a result, a gap occurred 
where a general purpose flight simulation program operating on low cost graphics 
workstations, is needed. A system that allows the flight characteristics of any aircraft be 
modeled on a low cost, general purpose, graphical workstation, would be of great use to 
aerodynamic students in studying the stability and control of different aircraft. 

B. FOCUS 

There are several issues to be addressed when incorporating an aerodynamic model 
into a computer simulation. The complexity of the aerodynamic model, which orientation 
model to use and how aircraft data should be represented in the system ate the critical issues 
and are center of focus in this work. 

I. Complexity of the Aerodynamic Model 

Of primary importance is that the complexity of the model fit the objective of the 
simulation. A complete aerodynamic model that includes fully articulated control surfaces 
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and airflow divergence patterns over the aircraft would seriously affect real time 
performance on any computer. Such models are usually computed in non-real time on super 
computers and are not appropriate for use on low cost graphics workstations. On the other 
hand, modeling the dynamics of an aircraft kinematically so that the aircraft’s velocity and 
orientation are a linear result of control input, does not reproduce the nuances of aircraft 
motion and response that a user of a flight simulation would expect. The aerodynamic 
model's complexity must provide as much realism as possible without reducing the frame 
rate below an acceptable level. Therefore, the focus in creating the aerodynamic model is 
to create an aerodynamic model of sufficient complexity that provides as much 
aerodynamic realism as possible without adversely effecting real-time performance. 

2. An Appropriate Orientation Model 

The choice of rotational model is fundamentally limited to either an Euler angle 
or quaternion approach. This has been the subject of heated debate among computer 
scientists as to which rotation model is the most appropriate [Goldiez and Lin,1991|. 
Basically, either model can be used to direct orientation, but, depending on the type of 
vehicle being modeled, one method has certain advantages over the other. The primary 
focus in defining an appropriate orientation model is in detennining which provides the 
right tool for the job [Shoemake,1985]. 

3. Defining Multiple Aircraft Aerodynamics Characteristics 

A general purpose flight simulator, such as is to be incorporated into NPSNET, 
should be capable of simulating an unlimited number of different aircraft types. Because 
each type of aircraft exhibits its own specific aerodynamics and handling characteristics, it 
is desirable to change these characteristics depending on which type of aircraft is to be 
piloted. One solution is to base the aerodynamic model on the aircraft’s stability 
coefficients, inertial coefficients, and airframe specifications, all of which ate generally 
available by consulting aerodynamic stability and control textbooks. Stability coefficients 
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provide a very accurate model of aircraft flight behavior. However, care must be taken 
when modeling some of the newer generation fighters. To improve maneuverability, these 
aircraft have been designed aerodynatnically unstable. Their' stability coefficients reflect 
this instability. 

C. SUMMARY OF CHAPTERS 

Chapter II covers those considerations for incorporating an aerodynamic model into 
NPSNET. A detailed aerodynamic model utilizing stability coefficients is presented. 
Chapter III investigates the difference between Euler angle and quaternion representations 
for defining rotations. The advantages and disadvantages of both methods are discussed 
and a model is presented for integration into NPSNET. Chapter IV provides 
implementation details of the aerodynamic and orientation model into the prototype flight 
simulator. Chapter V covers conclusions and lists requirements and suggestions for future 
work in this area. 




Figure 1.1: Flow of Data and Summary of Thesis Chapters 
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II. AERODYNAMIC MODEL 



A. INTRODUCTION 

In modeling die dynamic behavior of an aircraft within the framework of a computer 
based flight simulation, one usually begins by creating a mathematical model. This model 
represents the relationship between an aircraft and its interaction with the air mass in which 
it operates. Depending on the complexity of the model, additional forces, such as those 
associated with propulsion and environment, and complex systems, such as flight controls 
and weapons, are included [Rolfe,1986]. 

How detailed this model becomes is dependent on how it will be utilized. Generally, 
aerodynamic engineers will break the model down into components and study the results 
of detailed, but partial simulations. In complex flight simulators, such as those utilized to 
train pilots and aircrew, the model becomes more complex due to the increased number of 
systems that must be simulated and the response speed that is required for proper feedback 
to the pilot. As a result, the model loses fidelity, with the amount of detail limited by 
hardware and real-time constraints. 

Flight dynamics in NPSNET requires a mathematical model similar to that found in 
complex flight simulators. The framework for such a model is presented in Figure 2.1 . By 
adjusting the detail of those functions which surround the basic mathematical model, a 
simulation is built which can be made more or less detailed based on the speed of the 
workstation on which it is implemented. 
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Figure 2.1: Basic Aerodynamic Model 



B. COORDINATE SYSTEM 

Coordinate systems and the method in which they are described vary to a great extent 
on the application and the preference of the user. In aircraft simulations, coordinate systems 
fall into two broad classes, “body” coordinates and “earth” or “inertial” coordinates 
[Rolfe,1986]. Body coordinates have their origin based at the center of gravity and 
continually move with the aircraft. Inertial coordinates, on the other hand, are defined with 
respect to the earth and have their origin positioned at some suitable location such as the 
center of the simulated world. Other coordinate systems exist, based on parameters such as 
the flight path and angle of attack. However, the aerodynamic model presented in this paper 
is based on geometric body and world coordinate systems. In general, all aerodynamic 
forces, accelerations and velocities are calculated in the body coordinate system first, and 
then converted to the world coordinate system prior to updating an aircraft’s position and 
attitude. 

Figure 2.2 shows the generally accepted convention for labeling of the axes in the two 
coordinate systems found in most aerodynamic textbooks [Anderson,1989]. Body 
coordinates are defined with the origin at the center of gravity (CG), the x axis along the 
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fuselage pointing out the nose of the aircraft, the y axis along the wing-line pointing out the 
right wing, and the z axis pointing out the bottom of the plane. World coordinates are 
defined with the origin based at a fixed point on the ground, the x axis pointing North, the 
y axis pointing East, and the z axis pointing down. Because of its limited effect, the 
curvature of the earth is ignored. 
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Figure 2.2: World and Body Coordinate Systems 



C. DEFINITION OF TERMS 

In a dynamics model, velocities, accelerations and forces are described in both world 
and body coordinate systems. Without an explicit description of the variables used to 
describe the model, confusion can arise. Most tenns described in this paper refer to the 
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geometric body axes. However if a reference is made to the world coordinate system the 
subscript “w" is used. See Figure 2.3 for terms defined in world coordinates. 



X w , ^ w> 


Aircraft location in world coordinates (feet) 


U W! v w , W w 


Aircraft velocity in world coordinates (ft/sec) 



Figure 2.3: Terms Defined within the World Coordinate System 

Terms without the “w” subscript relate to body axes and include linear and angular 
velocities, accelerations, forces, moments, angle of attack and sideslip angle. Figure 2.4 
provides a description of these tenns. Note that the direction of angular accelerations and 
velocities and moment terms are defined using the right hand rule around their respective 
axis (Figure 2.5). 



u, v, w 


Linear velocity along X, Y, and Z body axes (ft/sec) 


p, Q,R 


Angular velocity along X,Y, and Z body axes (rad/sec) 


v t 


Resultant velocity Vector (U 2 + V 2 + W 2 ) 


V e 


Wind velocity across tail of aircraft 


U, V, w 


Linear acceleration (ft/sec 2 ) 


P, Q, R 


Angular acceleration (rad/sec 2 ) 


F x , F y , F z 


Forces acting on aircraft 


L, M, N 


Moments about the X, Y, and Z axes 


a 


Angle of attack[tan _l (W/U)l 


P 


Sideslip [tan' 1 (V/U)] 



Figure 2.4: Terms Defined within the Aircraft Body Coordinate System 
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The aircraft control surfaces sucli as elevator, ailerons, and rudder ate defined as a 
rotation in radians around their respective hinge points on the aircraft. When a control 
surface is flush with the aircraft, the angle of deflection is zero. See Figure 2.6 for a further 
description. 




Figure 2.5: Notation with Respect to Body Axes 



8e elevator deflection positive down (radians) 

a positive 8e produces a positive lift and a negative pitch moment 

8a aileron deflection positive left (radians) 

a positive 8a produces a negative roll moment 

8r positive nose left a positive 8r produces (radians) 
a positive sideforce and a negative yaw moment 



Figure 2.6: Terminology Defining Aircraft Controls 
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Within the aerodynamic model, the particular aircraft being modeled is characterized 
by certain dimensional characteristics. A description of these tenns is included in Figure 
2.7. 



s 


surface area of wing (ft 2 ) 


b 


wing span (ft) 


c 


cord length (ft) 


w 


weight (lbs) 


^xx 


roll inertia (slug-ft 2 ) 


Ivy 


pitch inertia (slug-ft 2 ) 


^zz 


yaw inertia (slug-ft 2 ) 


Ixz 


roll-yaw cross inertia (slug-ft 2 ) 



Figure 2.7: Aircraft Dimensional Specifications 



The model presented in this paper utilizes English Standard Units of measurement. 
Therefore, the following units are utilized: 



Distance: 


feet 


Force: 


pounds 


Mass: 


slugs 


Angles: 


radians 


Temperature: 


Kelvin 


Time: 


seconds 



In NPSNET distances are measured by meters. A conversion can be made into the 
metric system by either changing all the units of measurement to standard metric units prior 
to calculating any result, or converting the final result to meters prior to updating position 
using a simple conversion factor of 0.3048 meters/feet. In this paper the latter method is 
utilized. 
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D. MODEL DESCRIPTION 



As indicated in Figure 2.1, the mathematical model presented takes forces, control 
inputs and aircraft specifications as inputs, generating linear and angular velocities in 
aircraft body coordinates as outputs. Based on a classical representation of linear 
aerodynamics and utilizing the total force Equations (2.18 - 2.20), all the forces associated 
with lift and drag are calculated utilizing aerodynamic stability derivatives [Roskam,1979]. 
Stability derivatives, first used by Bryan over a half-century ago, assume that all 
aerodynamic forces and moments can be expressed as a function of the instantaneous value 
of the perturbation variables [Nelson.1989]. The perturbation variables are the 
instantaneous changes from the reference conditions of translational velocities, angular 
velocities, control deflections, and their derivatives. For example, the tenn 5X/5u is the 
stability derivative defining the change in X force with respect to the change in forward 
speed. This derivative can be expressed in terms of a non-dimensional coefficient C Xu as 
follows: 



Figure 2.8 lists the non-dimensional coefficients used in this model, these coefficients 
are generally broken down into three categories, lateral, longitudinal and control. The 
longitudinal coefficients represent forces effecting the longitudinal axes of the aircraft, 
while the lateral coefficients represent forces affecting the lateral axes of the aircraft. Non- 
dimensional coefficients, generated in actual aircraft testing, are available for most aircraft. 
By using these coefficients in combination with the dynamics equations, it is possible to 
build a general use flight simulator. 




( 2 . 1 ) 



where 




( 2 . 2 ) 



and Q is the dynamic pressure l/2pv^. 
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Longitudinal Coefficients 




C Lo 


Reference Lift at zero angle of attack 


C Do 


Reference Drag at zero angle of attack 


°La 


Lift curve slope 


C Da 


Drag curve slope 


C Mo 


Pitch moment 


C Ma 


Pitch moment due to angle of attack 


°LQ 


Lift due to pitch rate 


C MQ 


Pitch moment due to pitch rate 


C . 
La 


Lift due to angle of attack rate 


C 

Ma 


Pitch moment due to angle of attack rate 


Lateral Coefficients 


c vp 


Side force due to sideslip 


c L p 


Dihedral effect 


c LP 


Roll damping 


C LR 


Roll due to yaw rate 


C Np 


Weather cocking stability 


C NP 


Rudder adverse yaw 


C NR 


Yaw dampening 


Control Coefficients 


C L6e 


Lift due to elevator 


C D8c 


Drag due to elevator 


^M8e 


Pitch due to elevator 


C L5a 


Roll due to aileron 


c L5 , 


Roll due to rudder 


C N6a 


Yaw due to aileron 


C N6r 


Yaw due to rudder 



Figure 2.8: Aircraft Specification Notation 



1. Aerodynamic Forces and Moments 

Using the non-dimensional coefficients, lift, drag and sideforce are calculated as 

follows: 



L' = 



c c 


~(Vt + AVe) - 


2 - 


Ct + C T oc + C n Q — - — + C . be + C T c 8e 

Lo La v^2Vt La 2Vt L6c 


Vx 





P vrs (2.3) 



D = 



C Do + C Da a + 



(Vx + AVe) 
Vx 



.2-1 



pVx S 



(2.4) 
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S F - [^ypP + ^ V6i^ r J j 



(2.5) 



Once lift, drag and sideforce are calculated, these forces are translated into forces 
along the aircraft X, Y, and Z axes as shown in Equations 2.6 through 2.8. The temis f ax> 
Fay and Faz represent the resultant aerodynamic forces. 



r AZ = - V cos a- Dsina 


(2.6) 


: AX = L’sina- Dcosa-S F sinp 


(2.7) 


A Y = SpCOSp 


(2.8) 



The aerodynamic moments represent the torque forces about the center of the 
aircraft and are determined in the following equations 



L a = 



m a = 



n a = 



C LpP + C LP P 2Vx + C LR R 2Vt + + C LSr^ r 



C Mo + C Ma a+C MQQ2Vt +C M«^'2Vx +Clvl6e5e 



pV Sb 



(Vx + AVe) 
Vx 



C NpP + C NP P 2 Vx + CnrR 2VX + C N6a^ 3 + C N6r^ r 



pVx 2 Sb 



(2.9) 

p v x 2 Sc (2.10) 

2 

(2.11) 



The forces and moments that result from the above calculations are added to other 
forces and moments at this time (Equations 2.12 - 2.17). 



F X - F AX + F Thrusl 


(2.12) 


f y = f ay 


(2.13) 


II 

'"H 

> 

N 


(2.14) 


^ ^A F ^Torque 


(2.15) 


M = M a + M Thr|)st + M Gyro 


(2.16) 


N = N A + N Thrust + N Gyro 


(2.17) 
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Engine forces such thrust, torque anti gyroscopic effect as well as environmental 
forces such as wind shear can have anywhere from a minor to significant effect on the 
forces and moments along all axes of the aircraft [Roskam, 1979J. However, in order to 
limit the complexity of the model, some simplifications are made. Engine thrust is limited 
to the X-axis only and no calculations are made for torque or gyroscopic effect. This can 
easily be included at a future date if desired and when these values are available. 

2. Calculating Linear and Angular Acceleration 

The total force equations are used to determine the linear acceleration of the 
aircraft [Nelson, 1989J: 

F 

U = VR - WQ - gsinG + — (2.18) 

m 

P „ 

V = WP-UR + gsin<}>cose+ — (2.19) 

m 

p 

W = UQ - VP + gcos<j>cos9 + — (2.20) 

ni 

The total moment equations are used to derive the equations for solving for 



angular acceleration: 

L = IxxP- I xzR = I xz p Q+ (Izz-*yy)RQ (2.21) 

M = I yy Q + (l xx - I zz ) PR + I xz (P 2 - R 2 ) (2.22) 

N = ! zz p - ! xz p + (!yy - ! xx) p Q + ^zQR < 2 - 23 > 

However, prior to solving for either P or R, an interim step is required: 

L" = L + I xz p Q- Hzz-IyyJRQ (2.24) 

N' = N- (I yy -I xx )PQ-I xz RQ (2.25) 

Therefore, 

p = (^ Izz - ^ ^xz^ / (^xx^zz - ^xz^ (2.26) 
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(2.27) 



Q = <M- dxx-^zJPR-Jxz^-R 2 ))/^ 



R - ( N J xx + L *xz) ^ ^xx^zz J xz^ 



(2.28) 



are the equations for angular acceleration. 

3. Calculating Linear and Angular Velocity 

Linear and angular velocities are determined by numerically integrating the 
accelerations. The Trapezoidal Rule, sometimes referred to as the Modified Euler Method, 
is used. The general method for this integration technique is as follows: 



where 

P n = new value of P 
P n _, = previous value of P 

P n = current rate-of-change of P (obtained from Euler prediction) 

P n _ i = previous rate-of-change of P 
dt = integration step size 

4. Updating Position 

Because position updates occur in the world coordinate system, a system must 
exist to convert the aircraft’s linear velocities into changes in world position changes. 
Fortunately, this system was proposed years ago by Leonard Euler [EuIer,1758J. By 
knowing the current attitude of the aircraft in world coordinates, the linear body velocities 
are easily converted into world rates by applying the following transformation, which 
effectively rotates the [U V W] vector by the Euler angles[Roskum,1980]. The order of 
transformation is important when utilizing this method and proceeds as follows: 




(2.29) 



15 



v * 




10 0 
0 cos<f> -sin<J) 




U 

V 


w * 




0 sin<}> cos(J) 




w 



U, 

V 




cos0 0 sin8 




U *" 
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Integrating the resultant velocity vector, now in world coordinates, 
step of the program, a position update is obtained (Equation 2.33). 
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How Euler angles are obtained is addressed in the next chapter. 



( 2 . 30 ) 



( 2 . 31 ) 



( 2 . 32 ) 



by the time 



( 2 . 33 ) 
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III. METHODS OF ORIENTATION 



A. INTRODUCTION 

The aerodynamic model generates rotational velocities relative to the fixed aircraft 
body coordinate system. But, just as position updates could only be determined after 
converting linear velocities into world coordinates, orientation updates require conversion 
of angular velocities in a similar manner. Three methods exist for defining the conversion 
of angular velocities to orientations in world coordinates, each having its own particular 
advantages and disadvantages. 

The most popular of these three methods is known as the Euler Method. Using a 
sequence of three angles, the Euler Method provides an intuitive description of aircraft 
attitude in world space [Rolfe.1986]. These angles consist of the familiar azimuth angle y, 
the pitch angle 9, and the bank angle $. The next method, which has become popular in 
recent years, is the Quaternion Method. Based on the unit sphere, the Quaternion Method 
provides an elegant method of defining rotations through the use of four parameters. Three 
of the coordinates describe the axis of rotation while the fourth is determined by the angle 
through which the rotation occurs [Shoemake, 1985]. The third method of defining 
orientation is the Direction Cosine Matrix. Tire direction cosines relate the aircraft body 
axis frame to the world reference frame. Direction cosines, as used in flight simulation, can 
be detennined from either Euler angles or quaternions, are utilized for transformations 
between axes, and, alternatively can be directly updated by incremental rotation matrices 
[Paul. 1981]. 

Each method has its own particular advantages and disadvantages and their use 
depends on the application and the implementation. Because NPSNET is a networked 
simulator, the orientation model used must not only render orientations in the world of the 
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aircraft being piloted, but also of other aircraft in the world either flying autonomously or 
piloted remotely across the network. 

b. SURVEY OF METHODS 

1. Euler Method 

The most common method of defining an aircraft’s orientation in world space is 
by the Euler attitude angles. Starting with the aircraft’s axis origin aligned with the world’s 
axis origin, the Euler angles specify three successive rotations to bring the world 
coordinates into alignment with the aircraft. The fact that there exist 24 possible ways to 
define rotations, each with potentially different results, means that the order of these 
rotations is important. Euler chose the convention of rotating first about the z axis, then 
about the new x axis and finally about the newest z axis. This convention exists in celestial 
mechanics, applied mechanics, and molecular and solid state physics. The convention used 
in quantum mechanics, nuclear physics, and particle physics, chooses to rotate first about 
the z, then the new y, and finally the new z [Burchfiel,1990]. The convention most often 
used in graphics is standard to aerospace engineers and has been proposed for use by DIS 
[UCF/1ST 1990J. Using the right hand rule, rotations are made, first, counterclockwise 
about the z axis by the angle \|/, then counterclockwise about the new y axis by angle 0, and 
finally counterclockwise about the newest x axis by angle (Figure 3.1). 

The range of values the attitude angles can take are: 

n 

>l< = ±n d = ±~ ([> = ±jr 

The aerodynamic model generates velocities in body coordinates. As seen in 
Equations 2.30-2.33, the linear body rates were transformed into world rates by application 
of the Euler angles. What follows is a method for obtaining these angles. 
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There is a direct relationship between Euler attitude angles and the angular 
velocity of the aircraft around its body axes [Nelson, 1 ( 39], From this relationship 
(Equations 3. 1-3.3), the rates of change of the attitude angles can be derived. 



<j> = P + Qsin4>tan0 + Rcos<}>tan0 (3.1) 

0 = Qcos<t> - Rsim)) (3.2) 

\}' = Qsin<j)sec9 + Rcos<)>sec9 (3.3) 

The inverse of the above equations are 

P = <j> - \j/sin0 (3.4) 

Q = 9cos<)> + \j/sin<j>cos0 (3.5) 

R = — 0sin<{> + vgcos<{>cos0 (3.6) 



Equations 3.1 through 3.3 are also known as the gimbal equations and are quite 
commonly used in simulation. However, a problem exists when pitch, 9, goes through the 

vertical. That is, where pitch becomes +/-( k /2). At that point <j> and \j/ become undefined. 
Implementing a flight dynamics model capable of complete vertical maneuvering, 
necessitates “fixing” the code so a division by zero doesn’t occur. 
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2. Direction Cosines 



In the case of a flight simulation, transforming between body coordinates and 
world coordinates is done quite frequently. A convenient way to represent the 
transfonnation between two coordinate systems is with the direction cosines. Using matrix 
notation and the direction cosines (a, b, c), the transformation from body to world axes is 
expressed by: 



[ x w 




a, b, c, 
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Y w 
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^ 2 bj c 2 
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z w 




a 3 ^3 c 3 




Z_ 



(3.7) 



X, Y and Z represent vectors of any kind, such as force, velocity and acceleration. 
Hie inverse relationship, converting world coordinates to body coordinates is the 
transpose: 



X 




a 2 




x w 


Y 


= 


b | bo b 3 




Y W 


Z 




C 1 c 2 c 3 




z w 



(3.8) 



In terms of the Euler attitude angles, the direction cosines for the above 
transformations are shown in Figure (3.2) 



aj = cos0cos\|/ 


bj = sin <}> sin 0 cosy- cos<j>siny 


Cj = cos ()) sin 0 cosy ~ sin<J)siny 


= cos 0 sin y 


b^ = sin<}>sin0siny + cos<)>cosy 


C'y = cos<t>sin0sin y - sin^cos y 


= -sin0 


b 3 = sin(})cos0 


c 3 = cos<t>cos0 



Figure 3.2: Direction Cosines in Terms of Euler Angles 
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It should be noted that as there are 12 ways in which Euler angles can be defined, 
there also exist more ways to define the direction cosines. They depend on the coordinate 
system utilized. 

The direction cosines are absolutely necessary for transformations between 
coordinate systems, whether Euler angles or quaternions are used. Direction cosines w'ere 
already used for transforming linear velocities in body coordinates to world coordinates 
(Equations 2.30-2.32). By multiplying those transfonnation matrices into one matrix, the 
result would be identical to the transfonnation matrix of Equation 3.7. By using direction 
cosines, the need for detennining the intennediate velocities is eliminated. 

3. Quaternion Method 

An alternate method that has gained popularity in the graphics community in the 
mid 1980’s is through the use of the unit quaternion. Not a new' method, quaternions have 
been around for over a century. Augmenting the “Four-Parameter Method”, they have been 
useful to aerodynamic engineers for some time and are still the method of choice for 
describing spacecraft orientation [Mitchell, 1965]. Discovered by Sir William Rowan 
Hamilton in 1843 as a result of a search for a successor to complex numbers, quaternions 
provide an efficient means for updating orientations [Shoemake, 1985]. 

a. Defining Quaternions 

There are numerous ways to interpret the quaternion mathematically. The can 
be described as an algebraic quantity, 

w + ix + jy + kz (3,01 



as a point in three dimensional projective space (w, x, y, z), as a linear transformation of 
four space (Matrix), or as a scalar plus 3-vector: 



( w, v) v = ix + jy +kz 



(3.10) 
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The best notation depends on their intended use. The most intuitive approach 
is to view the quaternion as a scalar plus 3-vector (Equation 3.10). However, for algebraic 
manipulation. Equation 3.9, generally becomes more useful. 



A common way of defining quaternion orientation is in combination with 
Euler’s Theorem which states that the orientation of a rigid body can be described as a 
rotation about an axis v by rotation angle O (Figure 3.3)[Goldstien,1980]. Constraining the 
vector part to be of unit magnitude, the quaternion becomes: 

cos ^,v sin — (3.11) 

2 2 



This representation is always of unit magnitude such that: 

'y 1 ° 2 

w* + x + y + z = 1 



(3.12) 
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Prior to defining how to rotate a rigid body using the unit quaternion, it is 
necessary to review some of the mathematics associated with the quaternion. For the 
purpose of a flight simulation, an understanding of the multiplication and the quaternion 
derivative is necessary. More comprehensive reviews are obtainable from 
[Shoemake,1985; Golstein,1980]. 

b. Quaternion Multiplication 

Almost every application involving quaternions makes use of the 
mathematics associated with their multiplication. Similar to the algebra associated with 
imaginary numbers, quaternions have three imaginary units, i, j, and k and ate non- 
commutative under multiplication with 

i 2 = j 2 = k 2 = -1 (3.13) 



where 

ij = k = -ji jk = i = -kj ki = j = -ik (3.14) 

In algebraic notation, the product of quaternion Q multiplied by quaternion Q, is 



QQ, = (w + ix+jy + kz) (w, + ix, +jy, +kz,) 

= (ww, -xx, — yy,-zz,) 

+ i ( x\v j + wx j — zy j + yz j ) (3.15) 

+ j (yw, + zx, - vvy, + xz,) 

+ k ( zw , + yx , — xy | + wz , ) 



and is easily implemented in computer code. Vector notation is more intuitive to understand 
but not as simple to code: 

QQ| = (w,V)(w,,V|) = W W , - V • V | , W V | + W , V + V X V | (3.16) 

The result of the above multiplication is a rotation from the orientation represented by Q to 
the new cumulative orientation of Q and Q| in quaternion terms. 
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Multiplication provides a method of orientation extrapolation that can be of 
benefit in a networked simulation. If Q, represents a finite rotation based on an integral time 
step, and Q represents the cumulative rotation, then a repeated multiplying of Q by Q, will 
result in a smooth rotation across a series of update frames [Burchfield 990], 

c. Quaternions in Terms of Direction Cosines 

One frame of axes (body coordinates) can be brought into coincidence with a 
reference frame by a single rotation D about a fixed axis making angles A, B, and C with a 
second reference frame (world coordinates). The four parameters A, B, C, and D, therefore, 
define the orientation of the aircraft body in world coordinates [Rolfe, 19861. The 
transfonnation matrix relating the body to world coordinates using these four parameters is 
in Figure 3.4. 




By substituting the following quaternion parameters, this method is further simplified: 

q 0 = COS^D 

q, = cos A sin (3.17) 

q 2 = cosBsin~D 
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q 3 = cos C sin * D 

and the transfonnation matrix becomes, 

Xu, qo + qj- 42-^3 2(q,q 2 -q 0 q 3 ) 2(q () q 2 + q,q 3 ) 



W 

Vw 

Z w 



2(qiq 2 + q(iq3 ) q^i +c l: -( l3 2 (q 2 q 3 -q 0 q, > 
2 (qiq 3 -q 0 q 2 ) 2 (q 2 q 3 + q () q,) qy-q;-q 2 + q^ 







X 




Y 


2 


Z 


3 



(3.18) 



which represent the transform based on the unit quaternion. This result is used for 
converting determining position updates as well as orientation updates. 

d. Derivative of the Unit Quaternion 

To update the quaternion from angular acceleratHns, the following equations 

ate used: 

<io = (q,P + q 2 Q + q 3 R) 

qt = i(q 0 p + q 2 R-q 3 Q) (3.19) 

q: = ^(qoQ + q? p -qi R ) 

q 3 = \ (qo R + qiQ-q 2 p > 

Because of the constraint that the quaternion be of unit value and assuming an integration 
step size of less than 1, the above set of equations become: 

% = - l (qi p + q 2 Q + q 3 R ) + *.q 0 



qt = 2 (% p + q 2 R- q 3 Q) +M, 
q 2 = 2 (q«Q + q 3 p_ q i R ) + ^q 2 



(3.20) 



q 3 = 2 ( % R + q]Q-q 2 p ) +?t q 3 
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where 



* = 1 - <qo + < l| + l2 + ( ?3 ) 

e. Euler Angles from Quaternions 



(3.20) 



Many of the auxiliary computations involved with a flight simulation require 
the use of Euler angles. To obtain these angles from the transformation matrix (Figure 3.2), 
the following method is utilized. 

Because pitch is limited to ±n/2, the cos(0) is always positive. As a result, 
obtaining these angles is relatively simple. The elevation angle is derived as follows: 

Si " e = - ,a ’’ ,3.21) 

0 = asin ( a 3 ) 



To obtain the azimuth (y) angle it must be noted that, since the cos(9) is 
always positive, the sign value of a 2 always reflects the sign value of the sin(y): 

a, = cosGsiny (3.22) 

Therefore: 



“i 

V = acos (— g) • (sign [a 2 ] ) 
The roll (<J>) angle is similarly obtained: 

c 3 

<t> = acos( — ) • (sign [b 3 ] ) 



(3.23) 



(3.24) 



C. ADVANTAGES AND DISADVANTAGES 

Quaternions and Euler angles each have their own advantages and disadvantages. 
Euler tingles use only three components instead of four to represent orientation. If one were 
to send quaternions over a network in place of Euler angles, as has been proposed for DIS, 
network traffic would increase. However, in cases where angular rates remain constant for 
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long periods of time, by extrapolating orientation updates with quaternions, it would be 
necessary' to send an update only when an angle rate change occurs [Burchfiel,1990|. 

Quaternions can be computed directly from the dynamic equations, bypassing the 
computation of transcendental functions necessary in computing Euler angles. If each 
transcendental function costs approximately 20 arithmetic operations then the net cost for 
deriving a rotation update using the Euler Method is 94 operations [Burchfiel,1990]. This 
can be compared to 42 operations using the quaternion method, hi flight simulation, Euler 
angles are necessary’ for use in other simulator functions. This means that, approximately, 
another 64 operations are necessary, making the Euler Method slightly more 
computationally efficient. However, in many cases it is not necessary to compute the angles 
at each orientation update. Network updates should only occur when a rate change occurs. 
Additionally, radars, gyros and attitude indicators can be updated at slower rates. The 
bottom line is that orientations can be achieved very efficiently utilizing quaternions, 
without calculating Euler angles. When Euler angles are needed for other aircraft functions 
then quaternions become less efficient, depending on how often they need to be computed 
(Figure 3.5). 



OPERATION 


# EULER OPS 


# QUATERNION OPS 


DERIVATION 


96 


42 


CREATING ROTATION 
MATRIX 


20 


32 


TOTAL 


116 


74 


EULER ANGLE 
CONVERSION 


0 


64 


TOTAL CALCULATIONS 


116 


138 



Figure 3.5: Efficiency Comparison of Euler and Quaternion Methods 
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The most significant advantage of quaternions is that no singularity exists when pitch 
angle (0) passes through ±7i/2. In the Euler Method, P and R both become undefined due 
to division by zero. Techniques exist, however, for working around this singularity. 
Truncating values as ± 71 / 2 is approached will avoid this problem. If the pitch angle is 
truncated at values of 89.99 and 90.01 dien a 0.02 degree rotation skip results 
[Goldiez,1991]. Depending on the speed of the program and the rotation rates desired, this 
may not be noticeable. However, in more fine grained simulations, where a slow vertical 
maneuver is executed, it is a factor. 

D. PROPOSED ORIENTATION MODEL 

The use of quaternions is desirable for many reasons. Since a dynamics model is used 
to drive the piloted aircraft, and the graphics platfonn utilizes transfonnation matrices for 
rendering, using the quaternion method offers an efficient method of converting body rates 
to orientations in world space. 

Numerous aircraft operate autonomously within the prototype simulation that change 
very little in angular velocity. When this simulator is eventually networked to other 
workstations, quaternions will provide a way to foward interpolate the resulting changes 
in orientation on the remote computers, thereby eliminating the need for the continued 
transmission of update packets. If updates are eventually needed, the quaternion rates can 
quickly and easily be converted into Euler angles for transmission. 

As the number of aircraft increase in the simulated world, the number of calculations 
necessary for orientation updates begins to multiply. Since only the currently piloted 
aircraft makes use of the Euler angles for additional simulation functions, calculating Euler 
angles is not necessary for all other aircraft in the simulation. Therefore, it becomes more 
efficient to utilize quaternions for defining these rotations, saving approximately 42 
arithmetic operations per update per aircraft. 
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In the piloted aircraft, at the expense of 22 calculations per orientation update, a 
smooth rotation exhibiting no singularity about its axes during vertical maneuvers is 
possible. Future versions of this model will incorporate slow speed, fine grain functionality, 
so smooth changes in orientation about all axes is important and are best obtained using the 
quaternion. 
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IV. IMPLEMENTATION ISSUES 



Aside from converting the mathematical model directly to code, other considerations 
must be addressed in the simulation. First, it must be modular enough to be easily 
incorporated into the existing NPSNET networked battlefield simulation trainer. NPSNET 
requires that it be operable from a Silicon Graphics workstation using a spaceball 
peripheral for flight control input, and provide update information for transmission to other 
workstations utilizing existing network protocols. Second, it must be able to simulate many 
different types of aircraft and allow interaction with a potentially large number other 
aircraft, either piloted from other workstations or operating autonomously. Third, for future 
use by aeronautical engineers, the system must be designed so that the flight characteristics, 
of whichever aircraft being flown, are based on an aircraft’s dimensional characteristics 
and stability coefficients. 

A. OVERALL SYSTEM LAYOUT 



To satisfy the prototype’s basic requirements, the overall structure of the system 



is designed as shown in (Figure 4.1). An aircraft data file makes it simple to create new 




Figure 4.1: Prototype Flight Simulator Basic Structure 



aircraft, position these aircraft, and designate their handling characteristics. It is also be 
used to initialize the flight parameters of the autonomous aircraft. The flight data structure 
contains information describing the current state of the aircraft and its design 
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specifications. The aerodynamic model is used for updating the piloted aircraft’s body 
rates. The non-dynamic model also outputs a set of velocities, but unlike the dynamic 
model, it detennines these velocities in accordance with a predetennined script designated 
by the implemented. The orientation model converts body rates to position and orientation 
in world space. Euler angles are then detennined for the piloted aircraft and are needed to 
update cockpit displays. Autonomous aircraft do not require the calculation of Euler angles 
and bypass this function. 

B. THE AIRCRAFT DATA INPUT FILE 

Data records within the input file are divided into two categories, flight records 
and aircraft specification records. The flight record contains information describing the 
position of an aircraft and its initial flight parameters (Figure 1.2). The specification record 
contains the dimensional characteristics and stability coeff rents describing a particular 
aircraft. By manipulating the data in the specification record, one can change an aircraft's 
basic design as desired. To link an aircraft to a particular set of specification data, all that 
is necessary is to match the “type aircraft” information within the two records. This system 
allows one specification record to be used for an unlimited number of flight records. 

C. THE FLIGHT DATA STRUCTURE 

The flight data structure provides a global source of information on the state of each 
aircraft in the simulation. The infonnation contained here is necessary for operation of the 
aerodynamic and orientation models, and for updating cockpit displays (Figure 4.3). Note 
that the fourth item in this structure is a pointer to another data structure containing aircraft 
specification data. Not shown, this structure contains information identical to that found hi 
the specification record of the data input file (Figure 4.2). Maintaining flight information 
in two separate files saves some storage space by allowing more than one aircraft use a 
single set of specification data. 
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Flight Record 


1 


/id number/ 


1 


/type aircraft/ 


200.0 


/airspeed/ 


-930.0 


/posx/ 


300.0 


/altitude/ 


0.0 


/posz/ 


0.0 


/heading/ 


0.0 


/initial angle of bank/ 


3.0 


/initial gforce/ 


1 


/stalus:0-piloted,l-levelturn,2-g controlled turn/ 


Specification Record 


4 


/type aircraft-- A4/ 


1 


/jet or prop jet: 1 prop:0/ 


27.5 


M 


260.0 


/s/ 


10.8 


/c/ 


546.0 


'/nt/ 


8090.0 


Ax/ 


25900.0 


Ay/ 


29200.0 


M 


1300.0 


Axz/ 


0000.0 


/max thrust/ 


8000.0 


/mil thrust or horsepower/ 


0.03 0.3 


/CDo/CDa 


0.28 3.45 0.0 0.72 0.36 


/CLo /C La /CLq /CLda /CLde 


0.0-3.6-0.38 -1.1 -0.5 


/CMo /CMq /CMa /CMda /CMde 


-0.98 0.17 


/CYb /CYdr 


-0.12 -0.26 0.14 0.08 -0.105 


/CLb /CLp /CLr /CLda /CLdr 


0.25 0.022 -0.35 0.06 0.032 


/CNb /CNp /CNr /CNda /CNdr 


0.2618-0.5236 0.5236 


/deflection limits of nid, ail, elevator (radians)/ 



Figure 4.2: Data Input Records 
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Lnt 


ID: 


—number assignment 


int 


type: 


—aircraft type 


int 


status; 


—piloted:!) or autonomous level tum:l climbing tum:2 


typeptr 


Tptr; 


—pointer to aircraft specification data structure 


float 


Forces[3): 


—forces in X,Y,Zdir 


float 


Torques[3]; 


-torques around X,Y,Z axis 


float 


lincar_vel[3] : 


—velocity in X,Y,Z direction 


float 


angular_vel[3]; 


—angular velocity around X,Y,Z axes 


float 


linear_accel[3]: 


—linear accelerations 


float 


angular_accel[3]: 


-angular accelerations 


float 


sideslip: 


—sideslip or beta angle 


float 


ang_atk: 


—angle of attack 


float 


d_ang_atk; 


—angle of attack rate 


float 


lift: 


—total lift 


float 


drag: 


— total drag 


double 


Q14J: 


—quaternion 


Matrix 


H: 


—direction cosine matrix 


float 


euler_angles[31: 


-euler angles angles in radians - yawpitch.roll 


float 


pos[3]: 


—world position in X, Y, Z 


float 


ref[3J: 


-look direction 


float 


vel_world[3]; 


—velocities in world position - X,Y,Z 


float 


gfor: 


-amount of g force 


float 


rpm: 


—engine rpm 


float 


elev: 


-flight control positions 


float 


eltrini: 


-elevator trim 


float 


rud; 


—rudder position 


float 


ail: 


—aileron position 


float 


thro: 


—throttle position 


int 


flaps: 


—flap position 


int 


gear; 


—landing gear position 



Figure 4.3: Aircraft Flight Data Structure 



D. FLIGHT CONTROL INPUT 

1. Throttle 

The throttle in an aircraft is the pilot's primary means to control the engine. As 
such, a mapping of throttle position to engine rpm must be devised that incorporates delays 
associated with engine spool-up characteristics. In large, high speed simulators, rpm and 
engine data is retrieved from engine-specific lookup tables. Because tables such as these 
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are engine specific, a simpler, generic method was devised. By keeping track of the last 
throttle position, the following formula is utilized: 

Arpm = (T 0|d -T new )<pdt (4.1) 

where 

T = throttle position 
dt = delta time 

cp = engine spool-up delay factor 

Since applications using this simulator include both jet and propeller aircraft, it 
was decided that allowances should be made for the differing characteristics of the two 
types of engines. Jet engines are generally rated in terms of thrust (lbs), while propellers are 
rated in terms of horsepower (ft-lbs/s) [Anderson, 1989]. Since the aerodynamic model 
uses tlirust in tenns of lbs, it can use the data for jets directly. However, propeller driven 
aircraft require the following conversion: 

550t|HP p 

T = — — - — — (4.2) 

Vx p 0 

where 

r| = propeller efficiency (usually around .8) 

HP = engine rated horsepower 

= density altitude ratio where p o is the density at sea level 

2. Aileron, Elevator and Rudder Control 

A normal aircraft stick exhibits two degrees of freedom, left-right for aileron 
control and back-forward for elevator control. Therefore, control inputs from the spaceball 
were limited to these directions. The maximum deflection of the control stick is information 
entered via the specification records. It is a simple procedure to read deflection data from 
the spaceball and linearly map it to a control deflection somewhere between +/- max 
obtainable deflection. 
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Rudder input is handled in a similar manner to aileron and elevator control. 
Maximum rudder deflection is entered via the specification records so that its position can 
be linearly determined from the amount of rudder control input. Rudder input is via the 
keyboard using the left and right arrow keys. Although not the same as using rudder pedals, 
separating the rudder input from the control stick inputs forces the operator to be concerned 
with the issue of rudder and aileron coordination. 

3. CONTROL OF NON-PILOTED AIRCRAFT 



hi the current system, an autonomous aircraft can exhibit one of two behaviors, 
depending on the value of the “status” data in the flight data record. The first behavior is 
that of a level turn. The initial roll angle, also a part of the flight data record (degrees) 
provides the seed for calculating the desired body rates. From the bank angle, the amount 
of g force necessary to maintain level flight at this angle of bank is [Anderson, 1989]: 



G- force = 



cost 



(4.3) 



from the g-force, the rate of turn (ROT) in degrees/sec can be determined: 



ROT = 



gjG-force 2 ■ 



(4.4) 



where g is gravity (32.2 ft/sec 2 ). A level turn requires a balance between the angular 

velocities Q, pitch rate, and R, yaw rate. Therefore, from Equations 3. 4-3. 6: 

Q = ROT sintj) ^ 5) 

R = ROT cos<(> 



The second behavior is that of a non-level turn. This turn is based on the initial 
angle of bank and g force entered in the data input file. The yaw rate, R, is calculated the 
same as for a level turn. The pitch rate, Q, is also computed using these fonnulas but, 
instead of using the computed g force for a level turn, it is calculated using the amount of 
desired g force entered via the data file. The calculations necessary to compute the angular 
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velocities are completed during program initialization and do not impact on the run time of 
the simulation. 

E. SPEED OF AERODYNAMIC MODEL 

The aerodynamic model exhibits a tendency to “blow up" if the time step between 
calculations becomes too great. This becomes evident when the aircraft speed and rotation 
rates increase. The solution to this problem is to run the aerodynamic model at a faster rate 
than the rest of the system. The trick is to determine how fast. “Smart” integration schemes 
can be found in the literature that increase or decrease the number of time steps according 
to the amount of numerical divergence present in the integrated values [Press.et al,l990]. 

A simple method is used in this program to solve for this problem. Because there exists 
a direct relationship between an aircraft’s speed and its tendency to produce a numerical 
instability, the time step was adjusted in direct relationship with the aircraft’s airspeed. 
Prior to updating body rates in the aerodynamic model, aircraft airspeed is measured and a 
the number of time steps is calculated (Figure 4.4). The timestep factor used in the current 
program is 0.03 seconds. 

aerodynamic model 

read aircraft velocity 

time step for aero model = timestep_factor * velocity 
for computed_time_step loop 
do aero calculations 
update aircraft state variables 
end loop 

exit aerodynamic model 

Figure 4.4: Algorithm for Computing Time Step 



Increasing the time step does not decrease performance of the overall system. The 
speed of the aerodynamic model in the current implementation ranges from 17.6 ms to 18.8 
ms for time steps ranging from two to 160. Running at approximately 120 ms, the graphics 
pipeline definitely remains the limiting factor in this system. 
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V. CONCLUSIONS AND FUTURE WORK 



The tediniques presented in this paper have proven to be an effective method of 
implementing a graphical dynamic flight simulation on a matrix-based graphics computer 
in real-time. Like most research and academic projects, this aircraft simulator is structured 
to allow for the addition of more detailed functionality. Current work includes the 
integration of a weapons delivery system and avionics suite. 

The orientation model functions developed during the course of this paper will become 
pail of the C program library at the school, thereby providing an alternate and more flexible 
tool for manipulating solid objects in a graphical environment. Integration of this 
orientation model into other dynamic simulation systems is i o under investigation. 



37 



APPENDIX A 



/* FILE: aero.c */ 

/* AUTHOR: Jo© Cook© */ 

/* DATE: 1/9/92 */ 

/* DESCRIPTION: This fil© contains the mathematical model for all */ 

/* aerodynamic calculations. */ 

/* Input: aircraft state structure (P) */ 

/* aircraft specification stucture (T) */ 

/* Output : updated state structure (P) */ 

/***************************************************************************/ 



aero_model (acftptr P,typeptr T, float 

{ 

float Q, QS , QSc, QSb, c2u, b2u, g, w, m; 

float tot_vel; 

float densalt r at io , densalt ; 

float aoa, cos aoa, sin aoa; 

float cos_sideslip , sin_sideslip; 

float For_thrust_X; 

float ailpos , rudpos , elevpos; 

float Itemp, Ltemp, Ntemp; 

float du, dv, dw, dp, dq, dr; 

float diff, 

float loopdelta? 

float world [3 ] , vwor [3 ] 1 

float vaircraf t [3 ] , 

float hpower; 

float max__thrust; 

float loop iterations; 

int ixt, alt; 

Matrix A MATRIX; 



deltatime ) 

/♦misc. variables*/ 

/♦total velocity*/ 

/*atmosp©ric data*/ 

/*angle of attack data*/ 

/*sideslip data*/ 

/*thrust force*/ 

/^control positions*/ 

/*temporary terms*/ 

/*newly computed accels*/ 

/*temp var for rpm determination*/ 
/*deltatime / loop iterations*/ 
/*vel in world inertial frame*/ 
/*vel in aicraft frame*/ 
/♦horsepower*/ 

/*thrust */ 

/♦number of aero calc iterations*/ 
/♦index and loop vars*/ 

/♦temporary matrix*/ 



y****************^Qg|- for valid structure********************* / 
if (P == NULL | | T = = NULL) { 

printf ("****error: null pointer sent to ae ro . c* * * * \n\n" ) ; 

return; 



/****some initial assignments****/ 
m = T ->m ; 
g = GRAVITY; 
w = m * GRAVITY; 

P->gfor = P->lift/w; /*calculate current G's*/ 
if (P->u < 10.0) P->u = 10.0; 

tot_vel = fsqrt (P ->v* P -> v + P->w*P->w + P->u*P->u) ; 

/************ atmoshe ric density calculations******************/ 
if (P— >posy >= 0.0 && P— >posy < MAX_ALT) { 

alt = (int) (P->posy/ (ALT_INCREMENT*FT_TO_METERS) ) ; 
densalt * densalt array [ alt ] ; 
densaltratio ■= densalt /DENSITYSL; 

} 

else { 

densalt = DENSITYSL; 
densaltratio « densalt /DENSITYSL; 

} 



j *** ★■*•*****★** * *calculat© new rpm*****************************/ 
diff = (P->thro - P-> rpm) / (2 . 0/deltatime) ; 

P->rpm = P->rpm + diff; 

/*************** ca q cu j_ a t e new thrust**************************/ 
if (T->class == 0) { /* if a prop plane*/ 

hpower = PROP_EFFIC * T->horsepwr; 
max_thrust = (hpower * HPR2THRUST) / tot_vel; 
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For_thrust_X = max^thrust * P->rpm * densalt ratio; 

} 

else{ /*if a jet plane*/ 

if (T->max > T->mil) For_thrust_X = (T->max* (P->rpm) ) * densaltratio 

else For_thrust_X = (T->mil* (P->rpm) ) * densaltratio; 

) 



y********* *dete rmine control differentials********************/ 
ailpos = (P->ail * T->def_ail) ; 
rudpos = (P->rud * T->def_rud) ; 

elevpos = (P->elev * T->def_elev) + P->eltrim; 

/************ ac jjust time step for aero calculations********/ 
loop_iterations - (float) ( (int) (t ot_vel*T ENTH_OF_AIRSPEED_F ACTOR) ) ; 
loopdelta = deltatime * 1 . 0/loop_iterations; 

/********** *i OC p to reduce time step for dynamic model*****/ 
for(ixt = 0 ; ixt< ( (int ) loop_iterat ions ) ; ixt ++){ 
if (P->u < 10.0) P->u = 10.0; 

tot^vel = fsqrt (P->v*P->v + P->w*P->w + P->u*P->u) ; 

P->sideslip = fatan ( (P->v/P->u) ) ; 
sin_sideslip= (sin (P->sideslip) ) ; 
cos_sideslip= (cos (P— >sideslip) ) ; 

aoa = fatan ( (P— >w/P— >u) ) ; /* + (5 . 0 /FAP_T0_DEG) */ 

cos_aoa = (fcos (aoa) ) ; 
sin^aoa = (fsin (aoa) ) ; 

P->d_ang_atk = aoa - P->ang_atk; 

P->ang atk = aoa; 

Q = 0 . 5 *densalt *tot_vel*tot_vel ; 

QS = Q *T->S ; 

QSc = QS*T— >c; 

QSb = QS*T->b; 

c2u = T->c/ (2 . 0*tot_vel) ; 

b2u = T->b/ (2. 0*tot_vel) ; 

/********** calculating aerodynamic force 3 *********************/ 
P— >li ft = ( T->CLo + 

T— >CLa * P— >ang_atk + 

T— >CLq * P— >q * c2u + 

T— >CLdalpha * P->d_ang_atk * c2u + 

T->CLde * elevpos ) * QS; 

P->drag = (T->CDo + 

T— >CDa * P->ang_atk ) * QS; 

P — > side force « (T->CYb * P->sideslip + 

T->CYdr * rudpos ) * QS ; 

P->Fx = (P->lift * sin_aoa - 
P->drag * cos^aoa - 
P->sideforce * sin_sideslip) ; 

P->Fy = (P->sidef orce * cos_sideslip) ; 

P->Fz = (-P->lift * cos_aoa - P->drag * sin_aoa) ; 

/********** calculating aerodynamic moments* * ******************/ 
P->L = (T->CLb * P->sideslip + 

T->CLp * P->p * b2u + 

T— >CLr * P->r * b2u + 

T->CLda * ailpos + 

T->CLdr * rudpos ) * QSb; 

P— >M = (T->CMo + 

T->CMa * P->ang_atk + 

T— >CMq * P— >q * c2u + 

T->CMda * P->d_ang_atk * c2u + 

T->CMde * elevpos ) * QSc; 

P->N = (T->CNb * P— >sideslip + 

T— >CNp * P->p * b2u + 
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T->CNr * P->r * b2u + 

T->CNda * ailpos + 

T->CNdr * rudpos ) * QSb; 

/* + + + * ★ c j e t errn i ne total moments and forces*'********** 1 *'* 11 '****** 
P->Fx = P->Fx + For thrust X; 



! * * * * * ■{- @mpo r a ry cal cul a t ions********************************* 
Ltemp = P— >L + T->Ixz * P->p * P->q - (T->Iz - T->Iy) * P->r 
Ntemp = P->N - (T->Iy - T-> lx) * P->p * P->q - T->Ixz * P->r 
ltemp = (T— >Ix * T— >Iz - T->Ixz * T->Ixz) ; 



/****** dete rmine Linear and Angular Accelerations************ 
du = ( P — > v * P->r - P — >w * P— >q + P->Fx/m - g * P->sptch ) ; 

dv = (P->p * P— >w - P->r * P— >u + P->Fy/m + g * P->sroll * P- 

dw - (P->u * P— >q — P->v * P->p + P->Fz/m + g * P->croll * P- 



dp = (Ltemp * T->Iz + Ntemp * 
dq = (P— >M - ( ( T-> I x - T->I z ) 

- ( (P->p * P->p - P-2 

dr = (Ntemp * T->Ix + Ltemp * 

/ + ** dete rmine velocities with 
P->u = P->u + ( ( P->udot + du 

P->v = P — > v + ( ( P->vdot + dv 

P->w = P->w + ( ( P->wdot + dw 

P->p = P->p + ( ( P->pdot + dp 

P->q = P— >q + ( ( P— >qdot + dq 

P->r = P->r + ( ( P->rdot + dr 



T->Ixz) / ltemp; 

* P->p * P->r ) 

r * P— >r) * T— >Ixz ) ) / T->Iy; 

T->Ixz) / ltemp; 

trapizoidal integration****/ 
) / 2.0) *loopdelta; 

) / 2.0) *loopdelta; 

) /2 . 0) *loopdelta; 

) / 2.0) *loopdelt a; 

) /2 . 0 ) *loopdelta; 

) / 2 . 0 ) *loopdelta; 



^/ * * * * * * ***record last acceleration***********************/ 
P->udot = du; P->vdot = dv; P->wdot = dw; 

P->pdot = dp; P->qdot *= dq;P->rdot — dr; 

} / * **end loop***/ 

/***** ^incrementally rotate quaternion*********************/ 
increment_quat (P->p, P->q, P->r , P ->Q, deltatime) ; 

/*****update rotation matrix from quaternion********* ******/ 
r ot at i on_mat r i x_f rom_qu at (P— >Q, P->H) ; 

/***** ©xt ract euler angles from matrix*********************/ 
euler_angles_f rom_mat rix (&P->sroll, &P->croll, &P->sptch, &P->cptch 

&P->shdg, &P — >chdg, P— >h__mat rix) ; 
euler_ang__f rm_rot_mat rix ( &P->roll , &P->ptch, &P->hdg, P->H) ; 

P->pt ch = P->ptch * RAD_T 0_D EG; 

P— >roll = P-> roll * RAD_T0_DEG; 

P->hdg = P— >hdg * RAD_TO_DEG ; 

/ **** ca l cu late new world velocity**************************/ 
vaircraft [0] = P->u; 
vaircraft [1] = P->v; 
vaircraft [2] *= P->w; 
t ranspose_mat rix (P->H, A_MATRIX) ; 

post_mult_mat rix_by_vector (A_MATRIX, vaircraft , world) ; 

vwor [ 0 ) « vworld [ 0 ] * FT_TO_METERS ; 

vwor [ 1 ] = -vworld [2] * FT_TO_METERS ; 

vwor [2] « vworld [ 1 ) * FT_TO_METERS ; 

/ * * * * calculate new world position**************************/ 

P — >H [ 3 ] [ 0 ] *= P— >posx += vwor[0] * deltatime; 

P->H[3][1] * P->posy +«* vworfl] * deltatime; 

P->H[3][2] = P->posz += vwor[2] * deltatime; 

P->H [3 ] [3] = 1.0; 



} 



y******** 
P->refx = 
P->refz = 
P->refy = 



calculate world reference point 
P->posx + (P->cptch * P->chdg) 
P->posz + (P->cptch * P->shdg) 
P->posy + P->sptch; 



* * 



************* 



* * j 



* * * / 

************ j 

* P->q ; 

* P— >q ; 

*************/ 

>cptch) ; 
>cptch) ; 
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APPENDIX B 



/•A:**********************************'***'**'********************************/ 



/* FILE: quaternions . c */ 
/* AUTHOR: Joe Cooke */ 
/* DATE: 1/9/92 */ 
/* DESCRIPTION: This file contains the functions for the quaternion */ 
/* rotation and matrix operations */ 



/★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★ ^/ 



/***multiply quaternion Qr by Q1 and return the product in Qr***/ 
mult_quate rnion (Ql,Qr) 
double Ql[4],Qr[4]; 

{ 

double n factor, q; 

n_f act or= 1.0/ { <Qr[0] *Qr [0 ] ) + {Qr [1 ] *Qr [1 ] ) + (Qr [2] *Qr [2 ] ) + (Qr [3 ] *Qr [3 ] ) ) ; 

Qr[0]=n_factor*(Ql[0]*Qr[0] -Ql[l]*Qr[l] -Ql[2]*Qr[2] -Ql[3]*Qr[3]); 
Qr[l]=n_f actor* (Q1 [1] *Qr[0] +Ql[0]*Qr[l] +Ql[2]*Qr[3] - Ql [3 ] *Qr [2] ) ; 
Qr [2 ] =n__f actor * (Ql [2 ] *Qr [0 ] + Ql[0]*Qr[2] + Ql[3]*Qr[l] - Ql [1 ] *Qr [3 ] ) ; 
Qr [3 ] =n_factor* (Q1 [3 ] *Qr [0 ] + Ql[0]*Qr[3] + Ql[l]*Qr[2] - Ql [2] *Qr [1 ] ) ; 

} 



/★★★convert roll, pitch and yaw angles to a single quaternion*****/ 
/★★★useful for intitalizing a quaternion with initial orientation***/ 
euler^to^quaternion (roll, pitch, yaw, Q)^ 
float roll , pitch, yaw; 
double Q [ 4 ] ; 

{ 

double r,p,y,cosr,cosp,cosy, sinr, sinp, siny ; 

r = ({ (double) roll/2.0) /RAD_T0_DEG) ? 

p = (double) ( (pitch/2.0) /RAD_TO_DEG) ; 

y = (double) ( (yaw/2.0) /RAD_TO_DEG) ? 

cosr = cos (r) ; cosp = cos (p) ; cosy - cos (y) ; 

sinr = sin (r) ; sinp = sin (p) ; siny = sin(y); 

Q[0] = (cosr*cosp*cosy) + (sinr* sinp* siny) ; 

Q[l] = (sinr*cosp*cosy ) - (cosr* sinp* siny) ; 

Q [ 2 ] = (cosr*sinp*cosy ) + (sinr*cosp* siny) ; 

Q [3 ] = (— sinr* sinp* cosy ) + (cosr* cosp* siny) ; 

) 



/* * *increment a quaternion from incremental body angular rates******/ 
/★★★p=roll q=pit ch and r— yaw deltas******************************/ 
increment_quat (p, q, r , Q, deltatime) 
float p, q, r ; 
double Q [ 4 ] ; 
float deltatime; 

{ 

double factor, qql , qq2 , qq3 , qqO, dqO, dql , dq2 , dq3 , p2 , q2 , r2; 

qqo = Q [ 0 ] *Q [ 0 ] ; 
qql = Q [ 1 ] * Q [ 1 ] ; 
qq2 - Q [ 2 ] * Q [ 2 ] ; 
qq3 = Q [ 3 ] *Q [3 ] ; 

p2 = (double) (p*0 . 5*deltat ime) ; 
q2 = (double) (q*0 . 5*deltatime) ; 
r2 = (double) ( r*0 . 5*deltatime) ; 
factor = 1 . 0- (qq0 + qql + qq2 + qq3) ; 

dqO = — ( Q [ 1 ] * P 2 + Q[2]*q2 + Q[3]*r2) + Q[0]*factor; 
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dql - (Q [ 0 ] *p2 + Q [2 ] * r2 
dq2 = (Q[0]*q2 + Q[3]*p2 

dq3 = (Q [ 0 ] * r2 + Q[l] *q2 
Q [0] = Q [ 0 ] + dqO / 

Q[l] = Q[ 1] + dql; 

Q[2J = Q[2] + dq2; 

Q [3 ] = Q[3] + dq3 ; 

} 



- Q[3]*q2) + Q[l]*factor; 

- Q[l]*r2) + Q [2] * factor; 

- Q[2]*p2) + Q [3 ] *f actor ; 



/^**matrix for conversion from worldrates to body rates***/ 

rotation_matrix__from^quat (Q,M) 
double Q [ 4 ] ; 

Matrix M; 

( 

double f , wf , xf , yf , zf , xxf , yy f , zzf , xyf , xzf , yzf , wxf , wyf , wzf ; 
f - 2.0/ (Q[0]*Q[0]+Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3]) ; 



xf = f * Q [1 ] ; 
yf = f * Q[2] ; 
z f — f * Q[3] ; 
xxf = x f * Q [ 1 ] ? 



yy f = y f * Q [ 2 ] ; 

zzf = zf *Q [ 3 ] ; 
xyf = x f * Q [ 2 ] ; 
xzf = xf *Q [ 3 ] ; 



yzf = yf*Q[3] 
wxf = xf *Q [ 0 ] 
wyf = yf *Q [0] 
wzf = zf *Q [0] , 

M [ 0 ] [ 0 ] = (float) (1.0- (yyf + zzf) ) ; 

M [ 0 ] [ 1 ] = (float) (wzf+xyf) ; 

M[0] [2] = (float) (xzf-wyf) ; 

M [ 1 ] [0 ] = (float) (xyf-wzf) ; 

M [ 1 ) [ 1 ] - (float) (1 . 0- (xxf+zzf) ) ; 

M [ 1 ] [2 ] = (float) (yzf+wxf) ; 

M[2] [0] = (float) (xzf+wyf) ; 

M[2] [1] - (float) (yzf-wxf) ; 

M[2] [2] = (float) (1.0- (xxf +yyf ) ) ; 
M[3] [0] = 0.; M [3 ] [1] = 0.; M[3] [2] 
) 



= 0.; M [ 3 ] [3] 



1 . ; 



/***this function "extracts"' sine and cosine data from a cosine matrix 
represented in conventional form. *********************************/ 
euler_angles_f rom_matrix (sroll, croll, sptch, cptch, shdg, chdg, M) 
float *sroll, *croll, *sptch, *cptch, *shdg, *chdg; 

Matrix M; 



{ 



float sp, cp, sr , cr , sy , cy ; 
sp = — M [ 0 ] [2] ; 
cp = f sqrt (1 . 0- (sp*sp) ) ; 
if (cp *== 0 . 0) ( 
sy = 0.0; 
cy = 1.0; 
s r = M [ 1 ] [ 0 ] ; 
cr = M [ 1 ] [1] ; 

} 

else ( 



sy = M[0] [ 1 ] / ( cp ) ; 
cy = M [ 0 ] [ 0] / (cp) ; 
sr = M [ 1 ] [ 2 ] / ( cp ) ; 
cr = M [ 2 ] [ 2 ] / (cp) ; 



} 

* sptch = sp;*cptch = cp;*sroll 



) 



/*s inpitch* / 

/*cospitch* / 

/*if pitch at +/- 90deg*/ 
/*sinyaw */ 

/*cosyaw */ 

/*sinroll*/ 

/*cosroll* / 



/*sinyaw */ 

/*cosyaw */ 

/*sinroll */ 

/*cosroll */ 

sr;*croll = cr;*shdg = sy;*chdg = cy 
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/* * *xhe following function will return the roll, pitch and yaw from a standard 
transformation matrix. This routine should work for a transformation in a coordinate 
system from body coords where the x y and z axes move in pitch roll and yaw 
respectively, to world coordinates where X Y and Z point North, East and Down 
re spect ive ly ***■******'*★**'*★'*****'*******************'********'***'***'*'**★*'*'*/ 

euler_ang_f rm_rot_mat rix (roll, pitch, yaw, matrix} 
float * roll , *pitch, *yaw; 

Matrix matrix; 

{ 

float cospitch, sinpitch, roll_fact, yaw_fact; 

/**determine pitch angle in radians and cospitch**/ 
sinpitch = -mat rix [ 0 ] [ 2 ] ; 

cospitch «= fsqrt(1.0 - (sinpitch * sinpitch)}; 

*pitch = - (fasin (matrix (0 ] [2 ])) ; 

/**determine roll angle in radians**/ 
roll_fact = matrix [2] [ 2 ] /cospitch ; 
if (cospitch != 0.0} 

{ 

if (mat rix [ 1 ] [2 ] < 0.0) *roll = - (f acos (roll_f act) ) ; 

else 

*roll = f acos (roll_f act ) ; 

} 

else 

* roll = 0.0; 



/**determine yaw angle in radians**/ 
yaw_fact = matrix [0 ][ 0 ] /cospitch; 
if (cospitch != 0.0) 



if (matrix [ 1 ][ 0 ] < 0.0} *yaw = - (f acos (yaw_fact ) } ; 
else 

*yaw = facos (yaw_fact) ; 

} 

else 

*yaw = 0.0; 

/***this routine assumes matrix is in graphics convention form. The robotics convention 
requires an additional matrix transposition prior to utilizing this function*********/ 



e u 1 e r_t o_s g i_a x i s_c o n ve rt (M^matrix , W_mat rix} 
Matrix M_matrix, W_matrix; 

{ 



W_matrix [0 ] [0] 
W_mat r ix [ 1 ] [ 0 ] 
W_matrix[2] [0] 
W_matrix [0 ] [ 1 ] 
W_matrix [ 1 ] [ 1 ] 
W_matrix [2 ] [ 1 ] 
W_matrix [0 ] [2] 
W_matrix [1 ] [2] 
W_matrix[2] [2] 
W_matrix [3 ] (0] 
W_matrix [3 ] [1 ] 
W_mat r ix [ 3 ] [ 2 ] 
W_matrix [3 ] [3 ] 
} 



M_mat rix [ 0 ] [0] ; 
M_mat rix [ 1 ] (0 ] ; 
M_mat rix [ 2 ] [ 0 ] ; 
-M_matrix[0] [2] 
-M_mat rix [ 1 ] [ 2 j 
-M_matrix[2] [2] 
M_mat rix [ 0 ] [ 1 ] ; 
M_mat rix [ 1 ] [1] ; 
M_mat rix [ 2 ] [ 1 ] ; 
M_matrix[3] [0] ; 
M_mat rix (3 ] [1] ; 
M_mat rix [ 3 ] [2] ; 
1 . 0 ; 



/* **Returns a vecter result ***************************************************** / 
void post_mult_matrix by vecto r (Mat rix M matrix, float * V array, float * Result} 

{ 

Result [0] = V_array [ 0 ] *M__mat rix [ 0 ] [ 0 ] + 

V^array [ 1 ] *M_matrix [0 ] [ 1 ] + V_array [2 ] *M_matr ix [ 0 ] [ 2 ] ; 

Result [1] = V_array [ 0 ] *M^matrix [ 1 ] [ 0 ] + 

V_array [ 1 ] *M_matrix [1 ] [1 ] + V_array [ 2 ] *M_matrix [ 1 ] [2 ] ; 

Result[2] = V_array [ 0 ] *M__matrix [2 ] [0 ] + 

V_array [ 1 ] *M_mat rix [2 ] [1] + V array[2]*M matrix [2] [2] ; 

} 
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APPENDIX C 



A. INTEGRATION OF AERODYNAMICS INTO NPSNET 

To complete the integration of the aerodynamic model into NPSNET the following 
steps must be accomplished: 

•Modify the vehicle model structure to fully describe a dynamic aircraft. 

•Create a method to assign aircraft characteristics on program initialization. 

•Add cockpit display as the appropriate user interface. 

•Modify input devices for aircraft control (simulate aircraft control stick). 

•Add aerodynamic model for control of the piloted aircraft. 

•Add orientation model based on quaternions for aircraft rotation. 

B. NEW VEHICLE MODEL STRUCTURE 

The current NPSNET has one model for all types of vehicles. Tanks, helicopters, jets 
and jeeps all move under the same sets of rules, and only one basic data structure for all 
vehicles was previously created. Because air vehicles operate in a different dynamic 
environment than ground vehicles, it is necessary to create an additional data structure for 
their description. A logical way to add additional data for an aircraft vehicle is through the 
use of a union data structure. Imbedding a union data structure containing the new 
aerodynamic model within the original data structure allows the insertion of additional data 
needed for a more complex model, without having to re-engineer the entire program with 

a completely new data structure (Figure A3.1). 

...(attached to Old Data Structure For Vehicle Model) 
union{ 

DEFVEH defaultveh: 

AIRVEH airveh; 

GNDVEH gndveh: 

SHIPVEH shipveh: 

SIMNETVEH sinitype: 

) veh; 

Figure C.l: Union to Be Inserted into Current Veliicle Data Structure 
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Two data structures are added to handle the aerodynamic model. Figure A3. 2 shows 
the structure containing the vehicle state information and Figure A3. 3 shows the structure 
containing the vehicle specifications. To connect the two structures, the pointer, 
TYPEPTR, was created. 



typedefFLIGHTSPECS * TYPEPTR; 
typedef struct 
/ 




i 

int totalacft; 


/*track total number of aircraft in simul*/ 


TYPEPTR Tptr. 


/♦pointer to aircraft type information*/ 


float 


forc[3]: 


/♦lbs: Forces in X,Y,Z direction*/ 


float 


torq[3J: 


/*ft-lbs: Torques around X,Y,Z axis*/ 


float 


lin_vel[3]: 


/*ft/sec: linear vel(u,v,w) in X, Y,Z dir*/ 


float 


ang_vel[3J: 


/*rad/sec: ang vel(p,q,r) around X,Y,Z axes*/ 


float 


Iin_accel[3]: 


/*fps2: Linear accels in X,Y,Z direction*/ 


float 


ang_accel[3J; 


/*deg/sec2;angaccel(p,q,r) around X,Y,Z axes*/ 


float 


sideslip; 


/*sideslip or beta*/ 


float 


ang_atk; 


/*angle of attack*/ 


float 


d ang^atk: 


/*angle of attack rate*/ 


float 


lift; 


/*total lift*/ 


float 


drag; 


/*total drag*/ 


double 


quat[4]: 


/*quatemion*/ 


Matrix 


h_matrix; 


/♦direction cosine matrix of world position*/ 


float 


sroll; 


/*sine of roll*/ 


float 


croll; 


/*cosine of roll*/ 


float 


sptch; 


/*sine of pitch*/ 


float 


cptch; 


/*cosine of pitch*/ 


float 


shdg; 


/*sine of heading*/ 


float 


chdg; 


/*cosine of heading*/ 


float 


euler_angles[3]; 


/*euler angles in radians - roll, pilch, yaw*/ 


float 


vel_world[3]; 


/*meters/sec: velocity in SGI coords X,Y,Z */ 


float 


hdg; 


/*world heading*/ 


float 


gfor; 


/*g force exerted from back stick*/ 


float 


rpm; 


/♦fraction: 0.0 1.0 engine rpm*/ 


float 


elev; 


/*position: -1.0 to 1.0*/ 


float 


eltrim; 


/♦elevator trim*/ 


float 


rud; 


/♦position: -1.0 to 1.0*/ 


float 


ail; 


/♦position: -1.0 to 1.0*/ 


float 


thro; 


/♦position: 0.0 to 1.0*/ 


int 


flaps; 


/♦boolean: 0-flapsup 1-flapsdown*/ 


int 


gear; 


/♦boolean: 0-gearup 1-geardown*/ 


1 AIRVEH; 



Figure C.2: Data Stucture Containing Aircraft State Information 
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typedef struct 

f 

int 


type: 




int 


class; 


/* J :jet 0:prop*/ 


float 


span; 


/*wing span */ 


float 


sarea; 


pwing area */ 


float 


cord; 


/*mean aerodynamic cord*/ 


float 


mass; 


/*mass */ 


float 


lx; 


/*x axis inertial term*/ 


float 


Iy: 


/* y axis inertial term*/ 


float 


Iz; 


/* z axis inertial temi*/ 


float 


Ixz; 


/*cross inertial term*/ 


float 


horsepwr; 


/*used for props*/ 


float 


max; 


/*max afterburner thrust*/ 


float 


mil; 


/*max non afterburner thrust*/ 


float 


CDo; 


preference drag coefficient */ 


float 


CDa: 


/*W-force drag coefficient due to aoa */ 


float 


CLo; 


Preference lift coefficient */ 


float 


CLd alpha; 


Pdelta alpha */ 


float 


CLq; 


Plift due to pitch moment*/ 


float 


CLa: 


plifl due to angle of attack*/ 


float 


CLde; 


Plifl due to elevator deflection*/ 


float 


CMo: 


Ppitch moment coefficient at 0 alpha */ 


float 


CMq; 


Ppitch moment coefficient due to pitch rate */ 


float 


CM a: 


Ppitch moment coefficient due to aoa */ 


float 


CMda: 


Ppitch moment coefficient due to delta aoa */ 


float 


CMde; 


Ppitch moment coefficient from elevator motion*/ 
PY-force coefflcent due to side slip angle */ 


float 


CYb: 


float 


CYdr: 


PY-force coefflcent due to change in yaw rate */ 


float 


CLr; 


/*roll moment coefficient due to yaw rate */ 


float 


CLp: 


/*roll moment coefficient due to roll rate */ 


float 


CLda; 


proll moment coefficient due to change in ail */ 


float 


CLb; 


/*roll moment coefficient due to side slip */ 


float 


CLdr: 


/*roll moment coefficient due to change in yaw */ 


float 


CNda: 


/*aoa change effect on yaw moment */ 


float 


CNb: 


Psidesiip change effect on yaw */ 


float 


CNp; 


Proll rate effect on yaw moment */ 


float 


CNr; 


/*yaw rate effect on yaw moment */ 


float 


CNdr; 


Pyaw acceleration effect on yaw moment */ 


float 


def_rud; 


Prudder deflection +/- in radians */ 


float 


def ail; 


/* aileron deflection +/- in radians */ 


float def elev; 

} FLIGHTSPECS; 


Pelevator deflection +/- in radians */ 



Figure C.3: Aircraft Specification Data Structure 



C. INITIALIZATION OF THE AIRCRAFT DATA STRUCTURE 

To give the system user a convenient method of assigning aircraft specification data 
to each aircraft in NPSNET, a data input file (aero.dat), together with “read_acftdat.c”, is 
used. One simply enters the data for a particular type of aircraft in the data input file and 
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specifies, in the type aircraft field, a particular vehicle type. The vehicle types have been 
previously defined in NPSNET for distinguishing between types of vehicles. If an aircraft 
in NPSNET has a type number that does not match a type number entered in the data file 
then a default type is assigned by way of the “read_acftdat.c” program. This default type 
provides the “easiest” aircraft to pilot. An unlimited number of types can be entered 
through this data file. Users should insure that the number of these records are specified at 



the top of the file (Figure A3.4). 



H 


/label 




1 


/# of aircraft types 




t 


/record type/ 




30 


/match to type of aircraft/ 




1 


/jet or prop jet: 1 prop:0/ 




27.5 


/ b/ 




260.0 


IS/ 




10.8 


Id 




546.0 


/m/ 




8090.0 


fid 




25900.0 


fly/ 




29200.0 


M 




1300.0 


/ Ixz/ 




0000.0 


/max thrust/ 




8000.0 


/mil thrust/ 




0.03 0.3 




/CDo /CDa 


0.28 3.45 


0.0 0.72 0.36 


/CLo /C La /CLq /CLda /CLde 


0.0 -3.6 - 


0.38 -1.1 -0.5 


/C Mo /CMq /CMa /CMda /CMde 


-0.98 0.17 




/CYb /CYdr 


-0.12 -0.26 


0.14 0.08 -0.105 


/CLb /CLp /CLr /CLda /CLdr 


0.25 0.022 


-0.35 0.06 0.032 


/CNb /CNp /CNr /CNda /CNdr 


0.2618 -0.5236 0.5236 


/deflection of rud, ail + stab 



Figure C.4: Aircraft Data Record 



“Read_acftdat.c” contains the code for linking the aircraft structures together. It takes 
the entire vehicle array, determines which vehicles are aircraft, matches the type records in 
the data file with the specified vehicle types, initializes the two data structures (above), and 
places an updated vehicle structure back into NPSNET. 



D. COCKPIT DISPLAY 



In order to have an appropriate interface between the user and the aircraft simulation, 
a rudimentary cockpit was devised and inserted in the code at the end of the procedure 
“graphicsproc()” found in jeep.c. 
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The projection of the cockpit is ortho2 and the following code was added: 

prefposition(0,XMAXSCREEN,0,YMAXSCREEN 

winsetljeepwin): 

acft_cockpit(vehpos); 

For control of the radar sweep the following global variables were added: 



float 


deltatime; 


/* frames per second*/ 


float 


radar_pos; 


/*keeps track of radar look angle*/ 


float 


rdr_posx; 


/*radar sweep x axis position*/ 


int 


reverse_dir; 


/*radar scan direction*/ 



The following definitions were placed in the “aero.h” file: 

LEFTSCOPE 48.5 
RIGHTSCOPE 59.0 
TOPSCOPE 14.0 
BOTTOMSCOPE 3.5 
SCANRATE 72.0 
SCANWIDTH 120.0 

SCOPEW1DTH (RIGHTSCOPE-LEFTSCOPE) 

HALFSCANWIDTH (0.5*SCANWIDTH) 

SCOPERATE (SCANWIDTH/SCOPEWIDTH ) 

CENTERSCOPE (LEFTSCOPE+(SCOPEWIDTH/2.0)) 

COMPASSRADIUS 5.0 
TEN (COMPASSRADIUS-0.5) 

THIRTY (COMPASSRADIUS- 1.0) 

E. INPUT DEVICE MODIFICATION FOR AIRCRAFT CONTOL 



1. Spaceball 

Control input for an aircraft is limited to movement of the Spaceball along the 
forward-aft axis and the left-right axis. Therefore, the method in which the Spaceball data 
is read was modified. To accomplish this, the procedure, newsbvals() in jeep.c, was 
separated by an ‘if’ statement into two cases: (1) an aircraft and (2) not an aircraft. From 
the Spaceball data array, values [0] and [2] are converted into aileron and elevator position 
data. The range of values is limited to -1.0 and 1.0. 

2. Keyboard 

The keyboard has been modified to add rudder, elevator trim and thrust input. The 
“>" and “<“ keys control the rudder input. The rudder can be moved left and right in .05 
radian increments, and maintains a value between -1.0 and 1.0. 
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Elevator trim is controlled by the ‘e’/’E’ and ‘r'/'R’ keys. ‘e'/’E’ trims the 
elevator surface up and ‘r’/'R’ trims the elevator surface down. Thrust is controlled by the 
V/’SYa’/’A’.and ‘d’/’D’ keys, and has a value between 0.0 to 1.0. The ‘s’/’S’ key moves 
the throttle higher in 0.05 unit increments, the ‘d’/'D’ key moves the throttle higher in 0.3 
unit increments, and the ‘a’/’A’ key moves the throttle down in 0.05 unit increments. 

F. INTEGRATION OF THE AERODYNAMIC MODEL 

The aerodynamic model used can be found in Appendix A and is incorporated as a 
separate compilation file “aero.c”. Within “jeep.c” a test is made as to whether the vehicle 
being controlled is a ground or air vehicle. If it is an air vehicle, program control moves to 
procedure movetheaircraft() found in “jeepmot.c”, where the procedure call, aero_model(), 
is made. The aerodynamic model functions exactly as described in the prototype simulator. 
After returning from procedure aero_model(), position upd; s are calculated and all the 
variables used by NPSNET are assigned new values. 

One aspect of NPSNET is the ability to switch from one vehicle to another through the 
planar map display. Since the vehicle starts out with an initial heading and airspeed, it is 
important to “remember” this information when assuming control of the vehicle. This is 
accomplished by initializing the aircraft state structure (Figure C.2) with airspeed and 
attitude data within the procedure switchveh() in file “dogsncats.c”. It is important to 
remember at this point that the aerodynamic model will convert all zero airspeed values to 
10 ft/sec. If the desire is to have an immediately flyable plane at the vehicle switch, then 
the airspeed should be initialized to at least 100 ft/sec. 

G. INTEGRATION OF ORIENTATION MODEL 

The orientation model can be found in file aero.c in two separate places. The first place 
is at the end of procedure aero_model(), where it is used to determine pitch, roll and yaw 
data for the driven aircraft. The second place is in the same file but in a separate procedure, 
rotate_translate_acft(). The procedure rotate_translate_acft() can be used in future 
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NPSNET improvements when a quaternion rotation model is desired for updating the 
orientation and positioning of all non-driven vehicles in NPSNET (Figure C.5). 



rotate_translate_acft (ACFTPTR P, float dtime) /*f ramerate*/ 

{ 

Matrix A_MATRIX; 
float vel_wor [3] ; 

/******i ncremen tally rotate quaternion*********************/ 
increment^quat (P->ang_vel [0] , P->ang_vel [ 1 ] , P->ang_vel [ 2 ] , 

P->quat , dtime ) ; 

/***** U pdate rotation matrix from quaternion* ******** ******/ 
rot at ion_mat rix_f rom_quat (P->quat, P->h_matrix) ; 

/★**** Q vtract euler angles from matrix*********************/ 
euler_angles_f rom_mat rix (&P->sroll, &P->croll, 

&P->sptch, &P->cptch, 

&P->shdg, &P->chdg, P->h_matrix) ; 
P->euler_angles [ 1 ] = ( f asin ( (P->sptch) ) ) ; 

P->euler_angles [0] = ( f as in ( (P->s roll ) ) ) ; 

if (P->croll < 0.0) 

if (P->sroll < 0.0) P->euler_angles [0 ] = (-PI - P->euler_angles [0] ) ; 
else P->euler_angles [ 0 ] = (PI - P->euler_angles [0] ) / 

P->euler_angles [2 ] = ( f as in ( (P->shdg) ) ) ; 

if (P— >shdg < 0.0) 

{ 

if (P->chdg < 0.0) P->euler^_angles [2] = PI - P->euler_angles [2] ; 

else P->euler_angles [2] = TWOPI + P->euler_angles [2] ; 

) 

else 

if (P->chdg < 0.0) P->euler_angles [2] — PI - P->euler_angles [2] ; 

/ **** **** C al CU late new world velocity in sgi coords**********/ 
t ranspose_mat rix (P->h_mat rix, A_MATRIX) ; 

post_mult__mat rix__by_vector (A_MATRIX, P -> 1 i n_ve 1 , ve l_wo r ) ; 
P->vel_world[0] = vel_wor[l] * FT__TO_METERS ; 

P->vel_world[l] = -vel_wor[2] * FT_TO_METERS; 

P->vel_world[2] = -vel_wor[0] * FT_TO_METERS / 

/****cal CU late new world position is sgi coords*****************/ 
P->h_mat rix [3 ] [0] += P->vel_world [0 ] * dtime; 

P->h_mat rix [3 ] [1] += P->vel_world [1 ] * dtime; 

P->h_matrix [3 ] [2] += P->vel_world[2] * dtime; 

P->h_matrix [3 ] [3 ] = 1.0; 

world lookat point*****************/ 

P->lookatpt [0 ] * P->pos[0] + (cosptch * coshdg) ; 

P->lookatpt [1 ] = P — >pos [ 1 ] + sinptch; 

P->lookatpt [2 ] = P->pos[2] + (cosptch * sinhdg) ; 

) 



Figure C.5: Procedure rotate_translate_acft() 



50 



APPENDIX D 




Photo 1: Passing over the Airport 
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Photo 2: Approaching an Aircraft Turning Through North 




Photo 3: Target of Opportunity 
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Photo 4 Closing in On T^o Aircraft 
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